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Abstract 

The objective of this work is to compare a high-order solver with a low-order solver for performing 
Large-Eddy Simulations (LES) of a compressible mixing layer. The high-order method is the 
Wave-Resolving LES (WRLES) solver employing a Dispersion Relation Preserving (DRP) scheme. The 
low-order solver is the Wind-US code, which employs the second-order Roe Physical scheme. Both 
solvers are used to perform LES of the turbulent mixing between two supersonic streams at a convective 
Mach number of 0.46. The high-order and low-order methods are evaluated at two different levels of grid 
resolution. For a fine grid resolution, the low-order method produces a very similar solution to the high- 
order method. At this fine resolution the effects of numerical scheme, subgrid scale modeling, and 
filtering were found to be negligible. Both methods predict turbulent stresses that are in reasonable 
agreement with experimental data. However, when the grid resolution is coarsened, the difference 
between the two solvers becomes apparent. The low-order method deviates from experimental results 
when the resolution is no longer adequate. The high-order DRP solution shows minimal grid dependence. 
The effects of subgrid scale modeling and spatial filtering were found to be negligible at both resolutions. 
For the high-order solver on the fine mesh, a parametric study of the spanwise width was conducted to 
determine its effect on solution accuracy. An insufficient spanwise width was found to impose an 
artificial spanwise mode and limit the resolved spanwise modes. We estimate that the spanwise depth 
needs to be 2.5 times larger than the largest coherent structures to capture the largest spanwise mode and 
accurately predict turbulent mixing. 


Nomenclature 


a 

speed of sound 

b 

mixing layer thickness 

E 

total energy per unit mass 

f 

frequency 

H 

tunnel height 

k 

turbulent kinetic energy 

L 

domain length 

M c 

convective Mach number 

N 

number of grid points 

P 

pressure 

R 

heat flux 

Re 

Reynolds number 

T 

temperature 

t 

time 

U 

mean velocity 

A U = U 2 -U 1 

reference velocity 
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u* 

dimensionless mean velocity 

Ui 

Favre filtered velocity 

U 

mean streamwise velocity 

u'u' 

Reynolds normal stress 

u'v' 

Reynolds shear stress 

U rms U U 

rms streamwise velocity 

V rms V V 

rms transverse velocity 

w rms = 'Jw'w' 

rms spanwise velocity 

Xi 

Cartesian coordinates 

x, y,z 

Cartesian directions 

*ij 

spatially filtered shear stress 

t sgs 

i i) 

subgrid scale shear stress 

P 

spatially filtered density 

s* 

displacement thickness 

Pi 

laminar dynamic viscosity 
reference velocity 

Pt 

turbulent dynamic 

Superscripts: 


~ 

Favre-filtered 

- 

Spatially-filtered 

SGS 

subgrid scale 

Subscripts: 


t 

total condition 

1 

top stream 

2 

bottom stream 


I. Introduction 

Use of the Reynolds Averaged Navier-Stokes (RANS) equations has become standard practice for 
aerodynamic analysis. However the Reynolds averaging procedure introduces additional correlation terms 
that must be modeled. These models work reasonably well for the equilibrium wall-bounded and free 
shear layer flows for which they were calibrated, but often require retuning of model coefficients for other 
flows. These models also have severe limitations in predicting unsteady, separated, or transitioning flow. 
Increasing computational throughput has enabled the use of more sophisticated methods such as Direct 
Numerical Simulations (DNS) and Large-Eddy Simulations (LES). These tools are more deterministic 
and, unlike RANS, are less reliant on tunable constants. Consequently, they offer promising physical 
insights that RANS methods cannot. Despite this, there remains a great many issues to ensure the 
accuracy of LES. 

The primary objective of this work is to compare the benefits and shortcomings of a high-order 
method and a low-order method for performing LES of a compressible mixing layer. The high-order 
solver is the Wave-Resolving LES (WRLES) code (Refs. 1 to 3) that employs a Dispersion Relation 
Preserving (DRP) scheme. The low-order solver is the Wind-US code (Ref. 4) which employs the 
second-order Roe scheme. It is well known that second-order schemes are more convenient to use for 
complicated geometries due to their small stencil size. Additionally, second-order schemes can be easily 
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adapted to unstructured grids whereas high-order schemes frequently are limited to structured grids and 
simple geometries. The advantage of high-order methods is that they can achieve the same level of 
accuracy using significantly fewer grid points compared to the low-order schemes. The work herein sets 
to examine whether low-order solvers, given sufficiently fine enough resolution, can produce the same 
level of accuracy as high-order methods. This effort also attempts to address the current state-of-the-art 
practice in performing LES of a compressible mixing layer. 

Several technologies stand to benefit from an enhanced understanding of compressible mixing in free 
shear flows. Examples include supersonic combustion, exhaust nozzles, and internal flows present in jet 
engines and scramjets, to mention only a few. High speed turbulent flows are often investigated using 
Reynolds Averaged Navier-Stokes (RANS) models, and at times with dilitation-based compressibility 
corrections. It is widely known, however, that these compressibility corrections are empirical at best and 
do not replicate the actual physics of compressible turbulent mixing layers (Refs. 5 and 6). Therefore, the 
longer term objectives of this work are to better understand the fundamentals of compressible turbulent 
free shear flows using LES and to apply that insight towards improving existing RANS turbulence models 
for engineering applications. 

Brown & Roshko’s (Ref. 7) classical work on low-speed mixing layers made a key discovery in its 
noteworthy evidence of coherent large-scale spanwise oriented vortical structures. The formation of these 
large-scale “rollers” is triggered by the Kelvin-Helmholtz instability. As the spanwise rollers convect 
downstream they combine through pairing and tearing processes. This work led to a series of 
experimental studies investigating compressibility effects such as Chinzei et al. (Ref. 8), Papamoschou & 
Roshko (Ref. 9), Goebel and Dutton (Ref. 10), Samimy & Elliot (Ref. 11), Hall et al. (Ref. 12), and 
Clemens & Mungal (Ref. 13). The experimental studies observed that the most significant effect of 
compressibility on a turbulent mixing layer is suppression of the growth rate relative to an incompressible 
mixing layer at the same velocity and density ratios. The experiments demonstrate that compressibility 
effects are greatest between convective Mach numbers of 0.5 and 1. The particular case examined here is 
the turbulent mixing between two supersonic streams at a moderate convective Mach number of 0.46. 
Similarly, numerical simulations have been performed by several researchers including Foysi and Sarkar 
(Ref. 14), Freund et al. (Ref. 15), and Pantano and Sarkar (Ref. 16). The numerical simulations showed 
that dilatation terms are relatively small and not affected by increasing convective Mach number, whereas 
the pressure-strain terms in the Reynolds stress transport decrease monotonically with increasing 
convective Mach number. The redistribution of the Reynolds stresses in turn lowers turbulent production, 
which results in reduced mixing layer growth rate. 

In this work, LES is employed to numerically investigate a reference turbulent mixing layer at 
moderate convective Mach number, 0.46. Two different flow solvers are used. The low-order solver 
employs a second-order upwind scheme while the high-order method implements Bogey & Bailly’s 
(Refs. 17 and 18) 1 1 -point DRP scheme. The effects of numerical scheme, subgrid scale (SGS) modeling, 
and filtering are also investigated. Furthermore, a previous hybrid RANS-LES work concluded that an 
insufficient spanwise domain limited the development of the full three-dimensional structures (Ref. 19). 
Therefore, the work herein, in addition to its primary purpose of comparing high-order and low-order 
methods, also seeks to analyze the effects of increasing the spanwise extent while maintaining the 
spanwise resolution. This type of parametric study was recently used with the same low-order technique 
to investigate sensitivity of separated bluff body flows to spanwise width (Ref. 20). 

Ideally, LES would be performed not only in the mixing region but also along the upstream channels 
to resolve their turbulent boundary layers. Recently, Inoue & Pullin (Ref. 21) investigated the zero 
pressure gradient turbulent boundary layer through LES over momentum thickness Reynolds numbers of 
Re e — 0(1O 3 ) — 0(1O 12 ). This required a large number of grid points to resolve the boundary layer. In 
the present work focusing on the shear layer, the Reynolds numbers were Re e = 5.1xl0 6 and Re d = 

4.7 x 10 6 for the top and bottom channel flows, respectively. Fully resolving these turbulent boundary 
layers while also performing LES in the mixing region is computationally prohibitive. Therefore, the LES 
simulations are designed to focus on resolving the shear layer itself. This will also enable a side-by-side 
comparison between the high-order method WRLES with the low-order method Wind-US. 
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The paper is organized as follows: A description of the flow is provided in Section II. In Section III 
the formulation is discussed, while the details of the computational approach are given in section IV. A 
brief summary of the low-order solver Wind-US and the high-order solver WRLES is provided. Details of 
the grid, boundary conditions, and experimental conditions are also presented. Comparisons of the high- 
order and low-order method are presented in Section V. The section begins with a general description of 
the flow and then considers specifics of numerical method, SGS modeling, and filtering. The importance 
of spanwise width is highlighted. The results section finishes with a discussion comparing the low-order 
and high-order method again but at coarser grid resolutions. Conclusions are presented in Section VI. 

II. Description of the Flow 

The current simulations recreate case 2 of the experimental mixing layer tested by Goebel and Dutton 
(Ref. 10). The schematic in Figure 1 shows two supersonic streams separated by a thin splitter plate 
exiting into a long mixing section. The tunnel has a height of 48 mm and a width of 96 mm. The splitter 
tip separating the two flows has a thickness of 0.5 mm (an important length scale affecting vortex 
shedding) and the entrances of both streams into the mixing section have equal heights of 24 mm. The 
high-speed flow is located on the top with a Mach number of 1.91, while the bottom flow is at Mach 1.36. 
The convective Mach number calculated from the flow conditions is M c = 0.46, where M c = (U 1 — 
t/ 2 ) /(^i + a 2 )- The two flows are pressure matched at 49kPa. The total temperature of the top flow is 
578 K while the bottom flow is at 295 K. The velocities are U 1 = 700 m/s and U 2 = 399 m/s. 

The experiments used Laser Doppler Velocimetry (LDV) to measure velocities in the mixing section. 
The total viewing length of the optical instrumentation was 500 mm. Transverse velocity profiles at several 
streamwise locations were measured. The Reynolds normal stresses in the streamwise and transverse 
directions were measured; however, the normal stress in the spanwise direction was not measured. The LDV 
system was used to measure properties of the incoming boundary layers at the splitter tip. The boundary 
layer thickness, displacement thickness, and momentum thickness in each of the two streams feeding the 
mixing layer were measured just upstream of the splitter plate trailing edge. The spectra of the turbulent 
kinetic energy and the fundamental frequency of vortex shedding were not measured. 


^ Top Flow: M=l. 91, U= 700m/ s, P=49kPa 


48mm 


r\ 


500mm 


▼ Bottom Flow: M=1.36, U=399m/s, P=49kPa 


Figure 1 . — Schematic of splitter tip separating the top flow from the bottom flow. 


TABLE L— EXPERIMENTAL FLOW CONDITIONS 


Parameter 

Top flow 

Bottom flow 

M 

1.91 

1.36 

U (m/s) 

700 

399 

T t (K) 

578 

295 

P (kPa) 

49 

49 

P (kg/m 3 ) 

0.51 

0.79 
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III. Formulation 


The Favre-filtered Navier-Stokes equations are solved numerically with both the high- and low-order 
methods. The viscosity and thermal conductivity are allowed to vary with temperature. Favre filtering is a 
density-weighted spatial filtering procedure that is applied to the time dependent Navier-Stokes equations 
to remove small-scale fluctuations that are too small to be resolved. The continuity, momentum, and 
energy equations are: 


dp d , A 

aF + to- (p ' fi<) = 0 - 


( 1 ) 


d(pfti) 

dt 


d r -~ dp 

+ —( P u t u t ) + — 


d 


+ t?- gs 

IJ T L IJ 




( 2 ) 


^ (ujE + UjP ) = (uiiij + %f) + ( qj + qf GS ) • 


( 3 ) 


The density-weighted spatially-filtered velocity is denoted as u h while the spatially-filtered density and 
pressure are p and P, respectively. The stress tensor is f t j. For the scales that are too small for the scheme 
and grid to resolve, an optional subgrid scale model specifies rf GS and q GGS . 


IV. Computational Approach 

Two different numerical approaches are employed herein with the purpose of investigating whether a 
low-order scheme is capable of reproducing the results of a high-order scheme. High-order methods 
require less grid to resolve turbulent structures, whereas low-order methods require additional grid points 
to achieve the same level of accuracy. However, the large stencil sizes and complex numerics typically 
limit high-order methods to simple geometries. If the low-order method is capable of producing the same 
answer as the high-order method, then the low-order method may be more attractive for other geometries 
that are more complex than the relatively simple shear layer that is the focus of this paper. 

A. Flow Solvers 

1. WRLES 

The high-order solver is the Wave-Resolving LES (WRLES) code (Refs. 1 to 3). It is a finite 
difference code which employs the Dispersion Relation Preserving (DRP) scheme developed by Bogey & 
Bailly (Refs. 17 and 18). The stencil size is 1 1 points and is formally fourth-order accurate. While a 
traditional fourth-order scheme uses fewer points in the stencil, the additional points in the DRP scheme 
are used to minimize the dispersion error. It is an explicit central scheme similar to the classic tenth-order 
central scheme but with weights carefully chosen to minimize the dispersion error. Since it is a central 
differencing scheme, it is not inherently dissipative. A spatial filter, developed specifically to match the 
resolution of the differencing scheme, is used to preserve numerical stability. A coefficient multiplying 
the effect of the filter has been added, so that the magnitude of the dissipation can be controlled. The 
filter’s purpose and effects are strictly to help improve numerical stability. The DRP filter coefficient used 
was 0.5. Temporal discretization is performed using a low-dissipation and low-dispersion forth-order 
Runge-Kutta algorithm summarized in Berland et al. (Ref. 18). Three subgrid scale models have been 
implemented: Smagorinksy’s, Vreman’s, and a dynamic Smagorinsky. In particular, Vreman’s (Ref. 22) 
algebraic SGS model is used herein since it is less dissipative than the Smagorinsky subgrid model and 
not as computational intensive as the dynamic model. 
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WRLES uses the Message Passing Interface (MPI) standard to divide the domain across nodes and 
the Open Multi-Processing (OpenMP) directive to parallelize work on each node. In order to maintain 
order accuracy across interfaces, WRLES uses overlapping, point-matched zone coupling. 

2. Wind-US 

The low-order solver is Wind-US, the production flow solver of the NPARC Alliance. While it has 
historically been used primarily in RANS mode, the code may be used in LES or hybrid RANS-LES 
mode. It is capable of working with unstructured meshes as well as curvilinear structured meshes. It has 
several choices for temporal and spatial discretization as well as different options for flux splitting that are 
well documented (Ref. 4). For this particular problem, the numerical scheme used in the Wind-US CFD 
code is the second-order Roe Physical upwinding finite-volume scheme for spatial discretization and an 
implicit first-order time-stepping scheme. The Roe Physical scheme is based upon the original Roe 
scheme but modified for stretched grids. Note that while implicit temporal schemes require a matrix 
inversion, they are much more stable than explicit time schemes. When run in LES mode, Wind-US 
employs the Implicit (or Numerical) LES (ILES) philosophy. As discussed by Pope (Ref. 23), this 
approach assumes that the numerical scheme and grid act like a subgrid scale model, and so no explicit 
SGS model is used. 

For parallelization, Wind-US was run with the PVM (Parallel Virtual Machine) directive to create a 
parallel computing system. Wind-US uses abutting, point-matched zone coupling that is second-order 
accurate across interfaces. Within Wind-US, the MPI directive is also available for parallel computing. 

B. Grid 

The baseline grid used is structured as shown in Figure 2. The baseline mesh uses 14.4 million points 
to resolve the LES mixing section and has a thin spanwise width of 6 mm. Forty one points resolve the 
0.5 mm splitter tip yielding a Ay = Ax = 0.0125 mm in this region. Grid details such as the number of 
points, domain size, and minimum and maximum grid spacing are provided in Table 2 for the mixing 
region. The total length of the mixing section is L — 500 mm as may be observed in Figure 2. 
Downstream of the splitter plate towards the outflow the grid is stretched in the x-direction. The extent of 
grid stretching is deduced by comparing the minimum and maximum grid spacing. 



Inviscid Wall 



x 


T 

48mm 

_L 


500mm 


Figure 2. — Boundary conditions and dimensions of the LES mixing section (only every 8 th point is plotted). 


TABLE 2.— GRID PROPERTIES OF THE MIXING SECTION 


Parameter 

Direction 

X 

y 

z 

Number of points 

1025 

425 

33 

Length of domain (mm) 

500 

48 

6 

Minimum spacing (mm) 

0.0125 

0.0125 

0.1875 

Maximum spacing (mm) 

5 

0.2 

0.1875 
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C. Boundary Conditions 

1. General Boundary Treatment 

The flow conditions are taken from the experiments of Goebel & Dutton (Ref 10). The flow in the 
computational domain is supersonic throughout, thereby simplifying inflow and outflow boundary 
treatments (Figure 2). Inflows are fixed and the outflow is extrapolated. The upper and lower tunnel walls 
are assumed to be sufficiently far away from the shear layer such that they do not affect the mixing and 
are set as inviscid walls. Periodic boundaries are imposed on the spanwise sides of the domain. Their 
purpose is to simulate an infinite span and avoid the cost of modeling the sidewalls. The experiments had 
a finite span of 96 mm. As mentioned previously and discussed in greater detail later in this paper, the 
width of the spanwise dimension was investigated as a significant part of this work. 

2. Inflow Boundary Treatment 

In the experiment, turbulent boundary layers were present on the splitter plate upstream of the mixing 
section. Dutton & Goebel (Ref. 10) provided the values of the boundary layer thickness, displacement 
thickness, and momentum thickness. However, the actual velocity profiles and turbulent intensities were 
not reported. The best method for simulating these turbulent boundary layers would be a complete LES of 
the upstream region, producing unsteady turbulent boundary layers that match the experiment's integral 
thickness parameters. However, this type of analysis is cost prohibitive given the extremely small 
turbulent scales that would exist in a boundary layer at these Reynolds numbers. 

In this study, the inflow conditions were obtained from a separate two-dimensional RANS calculation 
of the mixing layer and upstream channels using the Menter Shear Stress Transport (SST) turbulence 
model. To do so, the lengths of the upstream ducts were varied until they produced RANS profiles that 
matched the experimental values of the displacement thicknesses. Table 3 compares the RANS-obtained 
profiles with the experiment by listing the computed values of displacement thickness 5*, mean 
streamwise velocity [/, and the area-averaged pressure P. Notice that the RANS solution indicates that a 
pressure-matched flow of 49 kPa is being delivered as in the experiments of Goebel-Dutton (Ref. 10). The 
incoming velocity profiles are depicted in Figure 3 and have displacement thicknesses matching the 


TABLE 3.— PROPERTIES OF THE INCOMING FLOW DELIVERED TO THE LES MIXING SECTION 


Solution 

Top flow 

Bottom flow | 

8* mm 

U,m/s 

P, kPa 

8* mm 

U,m/s 

P, kPa 

RANS upstream 

0.9318 

700.3 

48.898 

0.4681 

398.8 

48.940 

Experiment 

0.90 

700 

49 

0.44 

399 

49 



Figure 3. — Velocity profile of the flow entering the shear layer. 
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experiments to within 3.5 percent for the top flow and within 6 percent for the bottom flow. The resulting 
RANS mean-flow profiles were used as fixed conditions at the inflow plane to the mixing section. A 
limitation of this boundary treatment is that the RANS method only provides a steady state profile to an 
unsteady LES mixing region (Ref. 24). 

D. Time Averaging and Sampling Procedure 

To ensure numerical stability, the maximum time-step for the implicit low-order solver was found to 
be dt = 0.05 /is, while the explicit high-order solver was more stringent and the maximum time step was 
dt = 0.01ns. At the beginning of the simulations, there is a transient wave that travels through the 
domain. This originates from a normal shock produced at solution start-up. The time for all initial 
transient effects to be eliminated and before statistics can be gathered is no less than 5.5 Flow Thru Times 
(FTTs). The Flow Thru Time is defined as the time it takes a particle inside the shear layer to travel from 
the splitter tip to exit: 1 FTT = L/U AVG = 0.51 ms. U AVG is an estimate of the speed inside the shear 
layer taken to be the average velocity of the top and bottom streams. Extending the start-up wait time 
until 9.5 FTTs was found to be unnecessary, since statistics gathered after employing these two different 
start-up times were identical. 

The effect of sampling duration was examined by gathering second-order statistics over periods of 
9.5 and 33 FTTs that were then time averaged. Note that both of those sampling durations are 
considerably long. The resulting rms plots were compared and it was found that extending the sampling 
duration from 9.5 to 33 FTTs did not affect the rms values. Turbulent statistics were gathered at three 
frequencies corresponding to every 5, 1, and 0.4 ^5. Those frequencies were all extremely high and the 
results were found to be identical. A similar study varying sampling duration and frequency was 
performed by Mankbadi & Georgiadis (Ref. 20). After the turbulent statistics have been gathered, the 
power spectra is obtained by sampling the data every 0.05 /is after which a fast Fourier transform is 
performed to identify the fundamental vortex shedding frequency. 

V. Results and Discussion 

All the computations were performed on the NASA Pleiades supercomputing platform which is ranked 
amongst the top ten U.S. supercomputers. Pleiades is made up of tens of thousands of nodes capable of 
4.49 petaflop/s. Each node has two processors and each processor has either 6, 8, 10, or 12 cores. A typical 
simulation adhered to the following procedure: (1) initialize the flow field and wait for the start-up transient 
wave to exit the domain, (2) begin collecting data required for computing time-averaged turbulent 
intensities, and (3) at a much higher sampling rate (about one hundred times the fundamental) gather data to 
compute the power spectra of the flow. This procedure was used for both flow solvers. Wind-US 
simulations typically employed 39 nodes with a total of 468 cores and the turnaround time for the entire 
procedure was about a week. WRLES utilized 129 nodes with a total of 1536 cores and also had a 
turnaround time of about a week. 

Figure 4 shows dimensionless density contours of the mixing layer obtained using the two flow solvers 
for the baseline grid with 6 mm spanwise width. The high-order method, WRLES, is compared against the 
low-order method, Wind-US. Note that the size of the structures grows as the vortices pair and propagate 
downstream. There are more organized structures that can be seen in the WRLES plot whereas in the 
Wind-US plot the structures appear to be smeared (Figure 4(A) and (B)). The shock structure also slightly 
differs between the two. There are more reflections present in Wind-US than with WRLES. Wind-US’s 
small stencil and Riemann solver are adept at shock capturing whereas WRLES reduces the order to avoid 
smearing the shock with its larger stencil. The von Karman shedding is evident in Figure 4(C) and (D) 
where vortices are being produced just downstream of the splitter tip. The close-up views of the splitter tip 
demonstrate that the high-order method better resolves vortex shedding than the low-order method. Further 
downstream at x = 100 mm, the Wind-US solution appears to maintain the regular von Karman pattern 
while the WRLES solution appears to have broken down into a more irregular pattern. 
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0 5 10 x 15 20 250 5 10 x 15 20 25 

Figure 4. — Dimensionless density contours of the mixing section for (A) the high-order solver WRLES and 
(B) low-order solver Wind-US. Close up views of the splitter tip showing von Karman shedding for (C) WRLES 
and (D) Wind-US. 
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A. Flow Solver Comparison 

1. Effects of Flow Solvers 

Streamwise mean and rms velocities from the high-order method WRLES and the low-order method 
Wind-US are compared in Figure 5 (the transverse and spanwise rms velocities are presented in 
subsection B). Recall that the Wind-US approach employed a second-order finite volume solver and 
WRLES is a formally fourth-order accurate finite difference method employing an eleven-point DRP 
stencil. Both methods employ the ILES approach where no subgrid scale model is used. Despite the 
significant differences in numerical methodology, notice the relatively close agreement between the two 
numerical methods shown in Figure 5. The rms profiles at three different axial stations are very similar 
(Figure 5D to F), and the mean velocity profiles match fairly well (Figure 5(B) and (C)). The largest 
difference in mean velocity is in the wake at x = 50 mm in Figure 5(A). 

Figure 6 compares plots of the transverse velocity spectra of both numerical methods. The frequency 
is scaled by the fundamental frequency of vortex shedding f vs . The vortex shedding frequency, 
f vs , predicted is 124.51 and 134.28 kHz by Wind-US and WRLES, respectively. The data from the Wind- 
US solution shows a sub-harmonic exists at f/f vs = 0.5. The vortices produced at the tip recombine 
further downstream but at a lower frequency than that at which they were produced (Ref. 25). The 
vortices coalesce at the frequency of the sub-harmonic, f/f vs = 0.5. The sub-harmonic peak is not only 
present very close to the tip at x = 10mm but persists until x = 100mm. The fundamental frequency is 
the strongest at x = 10mm close to the tip and diminishes downstream, until it has disappeared by 
x = 100mm. 

The sub-harmonic cannot be seen in the spectra coming from WRLES. In addition, the spectra 
coming from WRLES shows that all the low frequency peaks have been damped out by x = 50mm and 
there is more energy in the smaller-scale, higher-frequency structures. The high-order numerics in 
WRLES can resolve smaller scales of turbulence, and these scales may be responsible for damping out 
the large-scale structures that persist in the Wind-US solution. Because of the differences in the 
fundamental frequency predicted and overall spectra, it is a little counterintuitive that two numerical 
methodologies would agree to such an extent in terms of the velocity and rms profiles shown in Figure 5. 

2. Effects of SGS Model 

In the previous section, the effects of the numerical scheme used by the high-order and low-order 
methods were examined. This section examines the effects of subgrid scale modeling on the 1 1 point 
DRP scheme solution. An additional case was computed with WRLES using Vreman’s subgrid model. 
The WRLES with SGS and WRLES without SGS cases are compared in Figure 7 and Figure 8. As 
evident from the data, the mean profile and rms profile are unaffected. The shedding frequency predicted 
by WRLES is 134.28 kHz for both cases and only minor differences are seen in the spectra. Therefore, 
the effect of SGS modeling is negligible for the 14.4 million point grid. In Figure 9, a plot of the SGS 
viscosity normalized by the laminar viscosity is provided. It is clear that the subgrid scale model is only 
active within the shear layer. Closer inspection of the scale reveals that /i t //i does not exceed unity 
anywhere in the mixing section. This explains why SGS effects are negligible for this problem. The grid 
is fine enough such that the vast majority of turbulence is computed directly and the SGS model does not 
add a significant amount of dissipation. If the grid was coarsened, then perhaps the SGS model would 
come into play. The fact that does not exceed unity also helps explain why the second-order method 
performed very well. The grid was fine enough in the mixing region, making any benefit of a higher order 
scheme irrelevant. The coarser grid and SGS effects will be explored in the next subsection. 
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Figure 5. — Comparison between the low-order solver Wind-US and the high-order solver WRLES on the 
14.4 million grid. For three different axial stations, the shear layer’s mean streamwise velocity (A) at 
x = 50mm, (B) at x = 100mm, and (C) at x = 150mm. The rms streamwise velocity (D) at x = 50mm, 
(E) at x = 100mm, and (F) at x = 150mm. 
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Figure 6. — Spectra of the transverse velocity for both 
numerical methods on the 14.4 million point grid. 
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Figure 7. — Effect of subgrid scale modeling on the shear layer’s streamwise mean and rms velocity at three 
different locations. The streamwise velocity (A) at x = 50 mm, (B) at x = 100 mm, and (C) at x = 150 mm. 
The rms streamwise velocity (D) at x = 50 mm, (E) atx = 100 mm, and (F) atx = 150 mm. 
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Figure 8. — Effects of SGS modeling on the spectra 
of the transverse velocity. 
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Figure 9. — The turbulent subgrid scale viscosity normalized by the laminar viscosity for the 14.4 million point grid. 
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3. Effect of Grid Density 

Recall in Section A-l, we demonstrated that the low-order solver could produce similar mean flow 
and turbulent statistics to the high-order solver on a fine grid. This section compares the two codes on a 
coarser grid. The coarsened grid used every other point in the streamwise and transverse directions while 
retaining the same spanwise width and points. 

The results of the grid study using the high-order solver WRLES are presented in Figure 10. Overall, 
the high-order solver performs nearly as well on the coarser grid. Close to the splitter tip, even the peak 
streamwise rms velocity values differ very little between the two resolutions. However, the mean flow 
predictions in the wake of the splitter tip are sensitive to resolution. Further away from the splitter tip, the 
mean flow and turbulent intensities predicted by the 3.6 million point grid are almost identical to the fine 
14.4 million point grid. Figure 10 also demonstrates that effects of SGS modeling are still very minor, 
despite the coarser grid employed. Comparing the SGS viscosities of Figure 9 for the finer grid to those of 
the coarsened grid in Figure 11, demonstrates that as the resolution was coarsened the SGS viscosity 
increases, but is still on the same order as the laminar viscosity. 

The results of the grid study using the low-order solver Wind-US are presented in Figure 12. These 
results demonstrate that the low-order method changes significantly at the coarse resolution. The low- 
order method cannot resolve the complex flow physics of vortex shedding, mixing, and shear layer 
growth. The inherent dissipation of the low-order scheme dampens out the instabilities and inhibits proper 
turbulent growth. There are substantial differences in both the mean flow velocities and the rms 
quantities. Not only is the peak rms value diminished but the extent of the profile across the shear layer is 
reduced significantly from that observed with the higher resolution. In other words, very little turbulent 
mixing is occurring. 

Figure 13 shows the transverse velocity spectra comparing the high-order and low-order methods at 
the coarse resolution. Note that the high-order method has well-defined sharp peaks at the fundamental 
frequency and its multiples. However, the low-order method has a blunt peak at the fundamental 
frequency. Recall the fundamental frequency predicted by WRFES on the fine grid was 134.28 kHz, here 
for the coarse grid it is higher at 141.91. Fikewise recall Wind-US predicted a fundamental frequency of 
124.51 kHz, for the coarser grid it is also higher 128.17 kHz. 

These findings demonstrate the expected result that high-order algorithms need less resolution than 
low-order methods to achieve the same level of accuracy. Unlike high-order methods, low-order methods 
are expected to fail at coarse resolutions and indeed did in this study. The increase in accuracy that high- 
order methods offer means that high Reynolds number problems become tractable at more reasonable grid 
sizes. With regards to computational wait time, high-order methods typically incorporate multi-stage time 
stepping and use large stencils that drives up the computational cost. Nevertheless, in many cases high- 
order methods can deliver that same level of accuracy at a cheaper computational cost since less grid is 
required. However, for the current problem, the two methods’ computational costs were about the same. 
This is primarily due to two reasons: (1) the explicit time-stepping of the high-order code requires a 
smaller time step than the implicit time stepping of the low-order method, and (2) the low-order method 
needs to regularly issue a checkpoint where the entire data is written to disk then an external process 
strips down the file to keep only the desired data. If not for this limitation in I/O performance, the low- 
order method would likely have been much faster than the higher-order method for the same grid. 
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Figure 10. — Grid study using WRLES with a spanwise width of 6mm. The mean streamwise velocity (A) at 
x = 50mm, (B) at x = 100mm, and (C) at x = 150mm. The rms streamwise velocity (D) at x = 50mm, 
(E) at x = 100mm, and (F) at x = 150mm. 
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Figure 12. — Grid study using Wind-US with a spanwise width of 6mm. The mean streamwise velocity 
(A) at x = 50mm, (B) at x = 100mm, and (C) at x = 150mm. The rms streamwise velocity (D) at 
x = 50mm, (E) at x = 100mm, and (F) at x = 150mm. 
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Figure 13. — Spectra of the transverse velocity for both 
numerical methods on the 3.6 million point grid. 


B. Effect of Spanwise Width 

Periodic boundary conditions are imposed in the spanwise direction to simulate an infinite span and 
avoid modeling the sidewalls. Selection of too small of a spanwise width cannot guarantee accurate 
capturing of the three-dimensional nature of vortex dynamics. To avoid imposing a nonphysical spanwise 
mode, the spanwise extent and resolution must be large enough to capture the wavelengths of the largest 
turbulent structures in the domain. 

An insufficient spanwise depth has been shown to artificially dampen the spanwise formation of 
vortices and resemble a two-dimensional instability (Ref. 20). In some cases, a two-dimensional 
instability predicts higher growth rates than a three-dimensional instability (Ref. 26). It is necessary to 
have a sufficient spanwise width to allow for the formation of three-dimensional coherent structures. Four 
simulations are performed here using the high-order method, with various spanwise widths in order to 
determine if: (1) artificial suppression is present when the spanwise width is smaller than the length scale 
of the largest coherent structure, (2) this artificial suppression is eliminated as the spanwise width is 
extended to be several times larger than the dominant length scale, and (3) as the mixing layer grows 
downstream, the size of the coherent structure grows thereby requiring a larger spanwise depth to be 
resolved. 

WRLES with the SGS model was used to study the effect of the spanwise width. Table 4 shows the 
grid resolution as the spanwise domain is doubled three times consecutively. Georgiadis et al. (Ref. 19) 
simulated this same flow and noted significant differences in rms statistics when examining much smaller 
spanwise depths ranging from 1 mm to 6 mm. The baseline grid was chosen to have the largest spanwise 
depth considered by Georgiadis et al. (Ref. 19) of 6 mm and had a spanwise grid spacing of A Z = 
6mm/32. Keeping the spanwise grid spacing fixed, the spanwise width was extended up to 48mm, 
utilizing a total of 1 12 million grid points. The peak values of w rms were also obtained. The streamwise 
and transverse resolution was unchanged. 


NAS A/TM— 20 1 5-2 1 874 1 


19 


TABLE 4.— SPANWISE WIDTH STUDY USING WRLES WITH SGS 


Spanwise width 

Resolution 

Grid points (million) 

Max Wrms /AU 

x = 50 mm 

100mm 

150mm 

6mm 

1025 x425 x 33 

14.4 

0.0947 

0.09343 

0.07341 

12mm 

1025 x425 x 65 

28.3 

0.1005 

0.1248 

0.1139 

24mm 

1025 x 425 x 129 

56.2 

0.1056 

0.1362 

0.1315 

48mm 

1025 x 425 x 257 

112.0 

0.1106 

0.1356 

0.1253 
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Figure 14. — Effect of spanwise width on the dimensionless streamwise vorticity. 

Figure 14 demonstrates the effects of increasing the spanwise width by plotting the streamwise 
vorticity at the three axial stations for spanwise widths of 6, 12, 24, and 48mm. Evidently the size of the 
structures is small at x = 50mm, but as they propagate downstream the structures grow in size and by 
x = 150mm the structures have grown considerably. The effect of increasing the spanwise width is most 
easily observed at x = 150mm. For a spanwise width of 6mm, the development of the structures is 
suppressed, but as the spanwise width is increased to 48mm, the structures are allowed to develop and 
thus appear to be bigger. The size of the turbulent structures relative to the mixing layer thickness can be 
seen by comparison with the value of b that is shown to the right of Figure 14. 

Although the mean flow is not affected by spanwise widths, the turbulence intensities are sensitive. 
The turbulence intensities in all three directions and turbulent kinetic energy are compared for all four 
spanwise widths at x = 50mm (Figure 15), x = 100mm (Figure 16), and x = 150mm (Figure 17). At 
axial stations of 50 and 100mm, the peak values of u rms and v rms are almost unchanged as the spanwise 
width was gradually increased from 6mm to 48mm. At x = 150mm (Figure 17(A) and (B)), the peak 
value in v rms is also unchanged but the peak u rms value slightly increases. From examining u rms and 
v rms , the need to extend the spanwise domain is not apparent since u rms and v rms are not changing 
appreciably. However, close inspection of w rms and k shows that w rms is artificially being suppressed for 
the 6mm case (Figure 15(C) and (D), Figure 16(C) and (D), Figure 17(C) and (D)). Expanding the width 
of the domain increases the peak w rms at the x = 50 mm and x = 100mm locations and k 
correspondingly increases. This differs from a similar study of spanwise width for a square cylinder 
computation (Ref. 20). For that problem w rms also increased with an increase in spanwise width but the k 
remained constant, forcing a rebalance of the energy between u rms , v rms , and w rms . 

Note that at x = 50mm, where the vortex size is small, the effect of the width on w rms is relatively 
small. However, as the turbulent structures increase in size moving further downstream, a larger spanwise 
width is required. Hence, at x = 100mm, note that the peak rms value is identical between the 24 and 
48mm case. Further downstream at x = 150mm, the peak value of w rms of the 48mm case lies between 
the 12 and 24mm cases. This demonstrates that extending the spanwise width eliminates artificial 
suppression of the spanwise rms value. 
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Figure 15. — At x = 50 mm, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 
(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 
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Figure 16. — At x = 100 mm, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 

(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 

In the aforecited work of Mankbadi & Georgiadis (Ref. 20) on a square cylinder, a spanwise depth ten 
times greater than the cylinder’s diameter demonstrated the w rms convergence. For a cylinder at low 
Reynolds number, the shed vortices are about the size of the cylinder’s diameter. However, for the shear 
layer examined here, the relevant length scale is the shear layer width, which grows downstream. The 
Kelvin-Helmholtz instability sets off the vortex shedding process that then becomes turbulent mixing 
downstream. As a consequence, the size of the largest structure increases axially. Hence, at x = 50mm 
where the length scale is small, a spanwise depth of 6mm is sufficient and spanwise width convergence is 
observed in Figure 15(C). Further downstream at x = 150mm, the size of the length scale has grown and 
hence a variation in w rms is observed in Figure 17(C). 

The relevant length scale is the shear layer thickness which correlates with the size of the largest 
coherent structure. A closer examination of the shear layer thickness is necessary to assess the ratio 
between the thickness and spanwise width at which spanwise width independence is reached. The shear 
layer thickness was defined by Goebel & Dutton (Ref. 10) as b = y 090 — y 0 10 , where y 090 is the 
location corresponding to U* = 0.90. Recall U* is normalized such that at the top flow it is one and at the 
bottom flow it is zero. From Table 5, note that at x = 50 mm the shear layer thickness is b — 2.5mm, 
while it has grown to the shear layer thickness of 4.4mm at x = 150mm (see Table 5). 
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Figure 17. — At x = 150 mm, the effect of spanwise width on the shear layer’s (A) streamwise rms velocity, 
(B) transverse rms velocity, (C) spanwise rms velocity, and (D) the turbulent kinetic energy. 


TABLE 5.— SHEAR LAYER THICKNESS 


Location 
x (mm) 

LES (z = 48mm) 
b(mm) 

Exp. 

b(mm) 

50 

2.549 

3.241 

100 

3.480 

4.508 

150 

4.426 

5.284 


Note that the u rms for the 6mm case in Figure 15(A) is identical to the 48mm case, thus spanwise 
width independence has been reached at 6mm. When the 6mm width is normalized by the 2.5mm 
thickness at x = 50 mm it yields a ratio of 2.3. Likewise in Figure 17(C) the peak w rms value for the 
48mm case lies between the 12mm and the 24mm cases, suggesting that they are relatively close at a 
spanwise width of 12mm. When the 12mm width is normalized by the 4.4mm thickness at x = 150mm it 
yields a ratio of 2.7. These results suggest that the solution becomes domain independent when the 
spanwise width is 2.5 times larger than the mixing layer thickness. 

For the shear layer case investigated herein, it was possible to extend the domain out until 48mm. 

This should ensure uninhibited growth of the turbulent structures, even at the furthest downstream 
measurement location of 450mm. For stations much beyond that, the spanwise width would likely need to 
be extended even more. However, at some point the spanwise domain cannot be extended without taking 
sidewall effects into consideration. The proposed factor of 2.5 cannot necessarily be generalized to other 
more complex flows. For other similar problems where resolution is scarce, a spanwise width that is only 
2.5 times larger than the relevant length scale may be sufficient to prevent artificial suppression of w rms . 
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VI. Conclusions 


This work compared a high-order method with a low-order method for performing large-eddy 
simulations (LES). The high-order method (WRLES) employed an eleven-point DRP scheme that is 
formally fourth-order accurate in space. For time stepping, it used the four-stage Runge-Kutta method, 
which is fifth-order accurate locally. For subgrid scale modeling, it uses Vreman’s model. The low-order 
method (Wind-US) is second-order accurate spatially and uses the Roe Physical scheme. For time 
stepping, it used an implicit backward Euler step that is first-order accurate locally, and does not utilize 
any subgrid scale modeling. EES were performed using both methods to capture the turbulent mixing 
between two supersonic streams at a convective Mach number of 0.46. The high-order and low-order 
methods were evaluated at two different levels of grid resolution. In order to compare with the 
experiment, the inflow profiles for LES of the shear layer were obtained using Reynolds Averaged 
Navier-Stokes (RANS) modeling over the upstream splitter plate such that the displacement thickness 
matched the experimental value reported. The results presented in this paper are summarized as follows: 

• For a fine grid resolution, the low-order method produced a similar solution to the high-order method. 
Both methods predict turbulent stresses that are in good agreement with experimental data. 

• Effects of subgrid scale modeling were found to be negligible at both resolutions using the higher 
order method and the Vreman SGS model. At even coarser resolutions, the SGS model might have a 
larger effect. 

• WRLES and WIND-US density contours were compared. For both methods, vortex shedding was 
observed to originate at the splitter tip and quickly transition to turbulence. More organized structures 
were seen with the high-order method whereas those structures appeared to be smeared with the low- 
order method. 

• The time for all initial transient effects to be eliminated and before statistics can be gathered is no less 
than 5.5 FTTs (Flow Through Times). The recommended sampling duration is 9.5 FTTs or longer. 
Any time-averaging window that meets these two conditions will not produce any substantial 
uncertainty. 

• When the DRP filter was nearly turned off, it was affirmed that its purpose is primarily for numerical 
stability and does not act as a SGS model. Thereby confirming that the reason the low-order and high- 
order method agreed was because the grid was sufficiently fine enough. 

• For a coarse grid of 3.6 million points, the high-order flow solver arrived at the same solution it 
reached with the 14.4 million point fine grid. However, the low-order solver did not demonstrate grid 
independence between the coarse and fine grid. Thus, high-order algorithms need less resolution than 
low-order methods to achieve the same level of accuracy. 

• We conducted a parametric study for the effect of the spanwise width on the accuracy of the solution 
by varying the width from a 6mm depth to a 48mm depth. As the shed vortices travel downstream and 
grow in size, the spanwise width may not be large enough to resolve all the spanwise modes and may 
introduce artificial ones. The peak root-mean-square of the spanwise fluctuation, w rms , was found to 
increase by fifty percent, when the spanwise width was increased from 6mm to 48mm. A spanwise 
width that is at least 2.5 times larger than mixing layer thickness was found to be necessary to capture 
three-dimensionality of vortex shedding. 

• Despite the differences in parallelization and time-stepping scheme, the two method’s computational 
costs were about the same. The larger stencil size and explicit time-stepping of the high-order code 
resulted in a larger computational cost for a given grid size compared to the low-order, implicit time- 
stepping code. However, the input/output performance of the low-order code was slowed by the need to 
regularly checkpoint the solution and spawn an external tool to extract the unsteady flow information. 
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The high-order method herein proved to be more accurate on a coarse mesh when compared against 
the low-order method. The small stencils of low-order methods are advantageous for their ability to 
handle complicated geometries whereas high-order methods are constrained by a large stencil size. The 
work herein establishes that low-order methods can be used for LES but require a larger number of grid 
points than their counterpart. Historically, most high-order methods do not retain their accuracy with 
unstructured grids. Recently, unstructured grids are favored for the inherent ability to redistribute mesh 
points where needed and remove unnecessary mesh points. Ideally the benefits of low-order and high- 
order methods would be combined by creating stable high-order unstructured algorithms. Such a method 
may push the state of the art by enabling high-order unstructured LES of practical applications. 
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